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The fixed node diffusion Monte Carlo (DMC) method has attracted interest in recent years as a 
way to calculate properties of solid materials with high accuracy. However, the framework for the 
calculation of properties such as total energies, atomization energies, and excited state energies is 
not yet fully established. Several outstanding questions remain as to the effect of pseudopotentials, 
the magnitude of the fixed node error, and the size of supercell finite size effects. Here, we consider in 
detail the semiconductors ZnSe and ZnO and carry out systematic studies to assess the magnitude 
of the energy differences arising from controlled and uncontrolled approximations in DMC. The 
former include time step errors and supercell finite size effects for ground and optically excited 
states, and the latter include pseudopotentials, the pseudopotential localization approximation, and 
the fixed node approximation. We find that for these compounds, the errors can be controlled to 
good precision using modern computational resources, and that quantum Monte Carlo calculations 
using Dirac-Fock pseudopotentials can offer good estimates of both cohesive energy and the gap of 
these systems. We do however observe differences in calculated optical gaps that arise when different 
pseudopotentials are used. 


I. INTRODUCTION 

Increasing demand for quantitative predictions of ma¬ 
terial properties, coupled with growing complexity of ma¬ 
terials of current research interest, promotes the devel¬ 
opment of high-accuracy methods from first principles. 
Recently the quantum Monte Carlo (QMC) method has 
emerged as a viable approach to quantitative calcula¬ 
tions of the properties of real materials with chemical 
identity^^®. QMC methods comprise a suite of tools 
for calculation of material properties via stochastic so¬ 
lution of the many-particle Schrodinger equation^*^^^^. 
Because of their direct treatment of electron correla¬ 
tion, QMC methods are amongst the most accurate avail¬ 
able and have a history of ground breaking, benchmark 
calculations^^. 

However, despite the successes of QMC thus far, an 
open question still remains as to the practical accuracy 
of the technique. In practice uncertainties arise from sev¬ 
eral approximations such as the use of pseudopotentials 
and the pseudopotential localization error^^ds^ finite size 
effects^^“^°, and the fixed node error that is present in 
fixed node diffusion Monte Carlo^^. At present, lim¬ 
itations in our understanding of the coupled effect of 
these uncertainties make practical usage of the technique 
challenging and hinder progress. The purpose of this 
manuscript is to present a systematic assessment of the 
uncertainties from these competing factors. As a case 
study, we have considered the semiconductors zinc se¬ 
lenide and zinc oxide in detail, and we quantify the size 
of uncertainties coming from the factors listed above. 

Variational Monte Carlo (VMC) and diffusion Monte 
Carlo (DMC) are two common approaches within the 
QMC suite of methods^°“^^. In variational Monte Carlo, 


an explicitly correlated form of the many-body wave 
function vL is used {e.g. Slater-Jastrow) and the ex¬ 
pectation value ('I'|iT|^')/('I'|'I') is evaluated stochasti¬ 
cally {H is the many-body Hamiltonian). The param¬ 
eters of the wave function can be optimized by either 
energy or variance minimization. In diffusion Monte 
Carlo (DMC), the wave function is described by a finite 
number of electron configurations (walkers). The time- 
dependent Schrodinger equation is mapped onto a dif¬ 
fusion equation in imaginary time, and the walkers are 
propagated according to the dynamics of the diffusion 
equation. This approach uses Green’s function methods 
to project out the ground state distribution of walkers 
from the trial wave function (with some complexities to 
be discussed later). Several recent demonstrations of the 
capability of QMC include the calculation of optical tran¬ 
sitions and thermal ionization levels for E-center defects 
in MgO^, intrinsic and extrinsic defects in zinc oxide^, 
metal to insulator transition in V02^, the volume col¬ 
lapse in cerium^, perovskite and post-perovskite MgSiOa 
in the earth’s lower mantle®, and several others®^®. 

Despite these promising studies, a clear set of “best 
practices” is not yet well-established for the for the calcu¬ 
lation of material properties such as total energies, band 
gaps, and defect properties within the fixed node DMC 
(DMC) framework. Understanding the relative size of 
uncertainties arising from competing factors like pseu¬ 
dopotentials, fixed node errors, and finite size effects will 
be an important aspect of making quantum Monte Carlo 
a standard computational tool. A recent comprehensive 
analysis of QMC applied to solids has demonstrated the 
viability of the technique, focusing on assessing the QMC 
prediction of lattice constants and equilibrium volumes 
across a extensive spectrum of materials (including met- 
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Approximation 

Controlled 

Description 

Assessment method 

Pseudopotential 

No 

Replace the core electrons with an ef¬ 
fective potential 

Vary within the space of reasonable 
potentials 

Localization 

No 

Approximate the diffusion Monte 
Carlo projector in the presence of a 
nonlocal potential 

Vary trial function and projector 
approximations 

Nodal 

No 

Fix the zeros of the trial wave func¬ 
tion when performing diffusion Monte 
Carlo 

Vary the trial wave function, apply 
variational theorem. 

Timestep 

Yes 

Approximate diffusion Monte Carlo 
projector 

Reduce time step until quantities are 
converged. 

Finite size 

Yes 

Finite size cells with periodic bound¬ 
ary conditions 

Increase the supercell size and average 
over twisted boundary conditions 


TABLE I. Sources of error in fixed node diffusion Monte Carlo that are considered in this work. The column “controlled” 
means that the error can be reduced to zero with a polynomially-scaling amount of computer time. 


als, semiconductors, and insulators)In this work, we 
take a complementary approach focusing instead on two 
materials in detail to assess competing sources of error 
in the calculation of total energies, atomization energies, 
and band gaps. While some aspects of this analysis may 
be known by experts, there are very few existing works 
that systematically assess these factors. This manuscript 
is intended to serve as a reference for those who wish to 
carry out these calculations. Also, even amongst practi¬ 
tioners there remains considerable discussion related to 
decoupling uncertainties arising from these effects. 

We chose zinc selenide and zinc oxide because they are 
both reasonably well understood, compound, wide band 
gap semiconductors from the II-VI family. The presence 
of zinc enables assessment of errors in the presence of 
localized states (the 10 3d electrons of Zn). The com¬ 
parison of the oxide to the selenide gives a chance to 
consider distinct effects arising from different chemistries 
since oxygen is more electronegative than selenium (ZnO 
is more ionic). As a semiconductor, ZnO has several 
desirable properties including good transparency, high 
electron mobility, and strong room temperature lumi¬ 
nescence. It is used as a transparent conducting elec¬ 
trode in liquid crystal displays and photovoltaic cells, and 
in electronics for thin-film transistors and light-emitting 
diodes. A large outstanding challenge is the elusiveness 
of obtaining p-type ZnO in the natively n-type material. 
For ZnSe, technological uses include II-VI light-emitting 
diodes (blue emission) infrared laser gain mediums 
(when doped with Cr)^^, infrared optical materials ex¬ 
hibiting a wide transmission wavelength range, and scin¬ 
tillators (when doped with Te)^®. 


II. METHOD 

A brief sketch of the quantum Monte Carlo methods 
used here follows. More details can be found in Refs [10- 
13]. We concentrate on explaining the methods well 
enough so that the approximations are clear; there are 
many details of implementation that affect the efficiency 


dramatically but do not affect the accuracy. Variational 
Monte Carlo is a direct implementation of the variational 
method for correlated wave functions. In this work, we 
use the Slater-Jastrow wave function, 

^'sj(R-) = Det[(()fe(rJ)]Det[^fc(r^)]exp[^ u(r*a,)], 

( 1 ) 

where (fik are one-particle orbitals obtained from a DFT 
calculation, (i,j) refer to electron indices, a is a nu- 
clear/ion index, R is the many-body electron coordinate 
(ri, r 2 ,..., r^r), and u is the same as the one in Ref. [26]. 
The parameters in the u function are optimized either 
using variance or energy minimization. The function u 
can take many forms, most common are two-body Jas- 
trow functions which explicitly include electron-electron 
interactions and three-body Jastrow functions which ex¬ 
plicitly include electron-electron-nucleus interactions. 

To obtain higher accuracy, we use fixed node diffusion 
Monte Carlo (DMC). In this method, starting with a trial 
function (in this work 'I't = 'I'sj), a projection to 
the ground state l$o) is performed: 

(R|$o)= lim / (R|exp[-iLT]lR')(R'l«'T)dR' (2) 

where H is the Hamiltonian. In principle, this integral 
can be done by Monte Carlo integration. In practice, 
there are two major impediments to the projection. First, 
the matrix element (R| exp[—i/r]jR') is only known in 
the limit of small r. This is solved by using the identity 

exp[—Fir] = I exp[—iFi5t]) and inserting many reso¬ 
lution of identity operators to increase the dimensional¬ 
ity of the integral. 6t is the time step in diffusion Monte 
Carlo. Second, the Hamiltonian matrix element includes 
a sign for fermions, which causes an exponentially de¬ 
creasing signal to noise ratio with system size. This is 
the sign problem, which can be resolved in general only 
with an approximation. We use the common choice of 
the fixed node approximation, in which $^jv(R') = 0 
wherever \I't(R-) = 0. The resulting method obtains 
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a variational upper bound to the ground state energy, 
with equality when the nodes of the trial wave function 
are equal to the nodes of the ground state wave function. 

While in principle DMC is a ground state method, the 
fixed node constraint allows approximate access to ex¬ 
cited states. Under some conditions, a variational upper 
bound to the excited state energy can be obtained^^, and 
in practice calculated excitation gaps can be quite accu¬ 
rate, even for challenging highly correlated materials^^. 
We perform excited state calculations by promoting an 
electron from the valence band maximum to the conduc¬ 
tion band minimum in the Slater determinant in Eq. (1) 
and proceeding as outlined for the ground state. 

In Table I, we present the major approximations 
present in the DMC technique, which we classify as ei¬ 
ther controlled or uncontrolled based on whether the un¬ 
certainty can be reduced to zero with a polynomially- 
scaling computer time or not. We will examine each of 
these one by one for the case of the semiconductors ZnSe 
and ZnO. In all cases for the work presented here, trial 
wave functions 4 't for the QMC simulations were gen¬ 
erated using the DFT framework as implemented within 
the CRYSTAL code^®. For the QMC calculations, we use 
the QWalk^® package. 

In our CRYSTAL calculations, we systematically opti¬ 
mize the localized basis sets using line optimization. The 
Gaussian exponents and weights are varied, and we se¬ 
lect the basis set that gives the minimum energy. The 
basis sets used here are available in supplemental mate¬ 
rial®*^. When this optimization scheme is used, our final 
DMC energies are not particularly sensitive to the basis 
set (the variation of the DMC energies are within error 
bars). Additionally, for all DMC simulations we use a 
sufficient number of walkers so that population bias is 
reduced to less than the stochastic error bars of our re¬ 
sults. In our simulations, we use the experimental lattice 
constants for the supercells of the materials in question. 


III. RESULTS AND DISCUSSION 
A. Uncontrolled Approximations 

1. Pseudopotentials 


TABLE II. Total energies (eV) of isolated atoms Zn, Se, and 
O according to DFT (PBEO) and DMC for both Trail-Needs 
(TN) and Burkatski-Filippi-Dolg (BED) pseudopotentials. 


Pseudo 

Atom 

DFT(PBEO) 

DMC 


Zn 

-1773.59 

-1744.48(3) 

TN 

Se 

-253.25 

-252.96(1) 


O 

-429.63 

-431.25(1) 


Zn 

-6179.75 

-6178.57(5) 

BFD 

Se 

-253.01 

-252.49(1) 


O 

-430.79 

-432.47(1) 


TABLE III. Atomization energy of zincblende ZnSe (top) and 
wurtzite ZnO (bottom) with experimental lattice constants 
according to DFT (PBEO) and DMC for both Trail-Needs 
(TN) and Burkatski-Filippi-Dolg (BFD) pseudopotentials, in 
comparison to experiment. 


Pseudo 

Method 

Atomization 
Energy (eV/fu) 

TN 

DFT(PBEO) 

5.71 

BFD 

DFT(PBEO) 

5.13 

TN 

DMC 

5.54(9) 

BFD 

DMC 

5.68(8) 

Experiment^'^ 

- 

5.511 


Pseudopotential 

Method 

Atomization 
Energy (eV/fu) 

TN 

DFT(PBEO) 

7.83 

BFD 

DFT(PBEO) 

8.63 

TN 

DMC 

7.61(8) 

BFD 

DMC 

7.67(8) 

Experiment“^ 

- 

7.59 


The optimal choice of pseudopotentials for quantum 
Monte Carlo calculations has emerged as an important 
question in recent years*^®’®^. In the pseudopotential ap¬ 
proximation, each electron is classified as either a core or 
valence electron, and the former are assumed largely in¬ 
ert, while the latter are significantly perturbed by bond¬ 
ing. The use of pseudopotentials eliminates the need 
to directly include the core electrons in the simulation 
and makes the calculation tractable, but it is an approx¬ 
imation and in reality a well-defined boundary between 
“core” and “valence” does not exist. 

We use relativistic Hartree-Fock {i.e. Dirac-Fock) 
pseudopotentials. There are several examples in the liter¬ 
ature that show they are well-suited for diffusion Monte 
Carlo simulations of solids®’^’®®’®®’®'*’. We consider Dirac- 
Fock pseudopotentials from two sources: the Burkatski- 
Filippi-Dolg®®’®® and Trail-Needs®'’’®® sets. These do not 
include core polarization or spin-orbit coupling. We focus 
our discussion on the zinc pseudopotential in particular, 
since the transition metal element with a full 3c?®® set 
of electrons is the more challenging and interesting case 
than O and Se by comparison. The Ne-core Burkatski- 
Filippi-Dolg Zn pseudopotential has 20 electrons in va¬ 
lence while the Ar-core Trail-Needs Zn pseudopotential 
leaves 12 electrons in valence. The comparison of the 
large and small core pseudopotential for Zn allows us to 
assess the extent to which the deeper semicore (Zn 3s 
and 3p) levels are perturbed by the bonding. There may 
be other differences between the two pseudopotentials as 
well, since they are constructed by different authors and 
not exactly in the same way. 

In Table II we show the DMC and DFT-PBEO total 
energies of the isolated Zn, Se, and O atoms for both the 
TN and BFD pseudopotential. As a test of pseudopo¬ 
tential accuracy, we show the atomization energy calcu- 
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lated within DMC in Table III for ZnSe and ZnO. For 
completeness we also show the atomization energies ac¬ 
cording to DFT-PBEO, but note that these do not neces¬ 
sarily give indications of pseudopotential accuracy, since 
the DFT prediction of the atomization energy itself may 
be wrong. The DMC total energies are obtained by fi¬ 
nite size supercell extrapolation and the T-moves scheme 
(both described in the next sections). The DFT results 
in Table III are obtained using the PBEO approxima¬ 
tion to the exchange correlation functional, and for the 
most part both the BFD and the TN pseudopotential give 
quite reasonable results within PBEO. For ZnSe, in com¬ 
parison to the experimental value of 5.51 eV/fu (formula 
unit)^^, the TN pseudopotential slightly overestimates at 
5.71 eV/fu while the BFD pseudopotential slightly under¬ 
estimates at 5.13 eV/fu. In both cases, DMC improves 
the description of the atomization energy: we obtained 
5.54(9) eV/fu for TN and 5.68(8) eV/fu for BFD. We ob¬ 
served similar trends in Table Ill for the atomization en¬ 
ergy of ZnO (7.59 eV/fu in experiment). Here, for both 
pseudopotentials DFT-PBEO results overbind the solid 
relative to isolated atoms: the degree of overbinding is 
0.24 eV/fu for TN pseudopotentials and quite large (1.04 
eV/fu) for BFD pseudopotentials. Once again, DMC im¬ 
proves the description substantially for both cases: TN 
gives 7.61(8) eV/fu and BFD gives 7.67(8) eV/fu. In all 
cases, the DMC results yield atomization energies that 
are within 0.1 eV/atom of the experimental value. 

For these materials, it appears that the large-core TN 
pseudopotential obtains similar results as the small-core 
BFD pseudopotential when using QMC, although we do 
not a priori expect this to hold true in general. Note that 
from Table Ill, it is not true for DFT (PBEO). As we will 
discuss later, it is also not the case for the DMC calcu¬ 
lation of the optically excited state. Atomic cores are 
more perturbed for highly ionic semiconductors, which 
may necessitate use of the small core pseudopotentials. 
Since ZnO is more ionic, for the remainder of the article 
we assess both TN and BFD potentials for that material, 
and only the TN potential for ZnSe. 


2. Localization error 

Nonlocal pseudopotentials introduce extra terms in 
the imaginary-time Green’s function (i?| exp[—17r]|i?') 
in DMC. These terms would cause a sign problem, 
which can be removed either using the localization 
approximation^® or the T-moves scheme^®. Both of these 
approximations introduce an error in addition to the 
fixed node error that depends on the quality of the trial 
wave function. When the trial wave function approaches 
the exact one, the localization error approaches zero. 

To assess the degree of the localization error, we gen¬ 
erated different trial wave functions within VMC by con¬ 
sidering (i) two and three body Jastrow functional forms, 
and (ii) energy and variance optimization of the Jastrow 
parameters. There are four possible combinations, and 


ZnSe, TN pseudopotential 



FIG. 1. Relative energy per formula unit ZnSe according to 
DMC, obtained by varying Jastrow factor and optimization 
method, both within the localization approximation and with 
T-moves. In the localization approximation, total energy dif¬ 
ferences can vary by as much as 0.75 eV for different Jastrow 
forms and optimization methods. By contrast, however, the 
variation is within error bars amongst all combinations when 
T-moves is used. The use of T-moves reduces localization 
errors, without substantial cost in the calculation time. En¬ 
ergies are shown relative to the minimum calculated DMC 
energy when T-moves is used. 


for all four trial wave functions we evaluated the total 
energy in DMC within the localization approximation 
and with T-moves. Fig. 1 shows the DMC energy per 
fu ZnSe using the TN pseudopotential for all four trial 
wave functions. The energies computed within the local¬ 
ization approximation (red bars) can vary by up to 0.75 
eV/fu from one scheme to the other, which shows a large 
dependence of the projected out wave function on the 
different forms of the trial wave function. By contrast, 
the incorporation of T-moves results in more uniform to¬ 
tal energies, which now vary within error bars instead 
(blue bars). Additionally, with T-moves the recovery of 
the variational theorem would allow unambiguous deter¬ 
mination of the best choice for the trial wave function, 
although in this case all four trial wave functions give 
statistically equivalent results. It is notable that the use 
of T-moves appears to reduce the dependence of the final 
DMC total energy on the trial wave function. 

In Fig. 2, we show the DMC relative energy per fu for 
ZnO for all four trial wave functions, but also compare 
the TN (Fig. 2a) and the BFD (Fig. 2b) pseudopoten¬ 
tial. The results for the TN pseudopotential in ZnO are 
similar to those of Fig. 1 for ZnSe: within the local¬ 
ization approximation the DMC total energies can vary 
by around 0.75 eV/fu whereas they are much more uni¬ 
form when T-moves is used. Interestingly, by contrast, 
in Fig. 2b for the small core BFD pseudopotential, the 
sensitivity of the total energies on the trial wave function 
is much smaller (in fact, within statistical error) within 
the localization approximation. In this case both the 
localization approximation and T-moves appear to give 
similar results, which is consistent with our expectation 
that localization errors are smaller when the pseudopo- 
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(b) ZnO, BFD pseudopotential 
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FIG. 2. Relative energy per formula unit of ZnO according to 
DMC, obtained by varying Jastrow factor and optimization 
method, both within the localization approximation and with 
T-moves. Within the localization approximation, localization 
errors are reduced with the small core BFD pseudopotential 
compared to the large core TN pseudopotential. Energies are 
shown relative to the miirimum calculated DMC energy when 
T-moves is used. 


tential core is smaller. 

For the remainder of the paper, all results presented 
have been obtained using T-moves and energy optimiza¬ 
tion to minimize the localization error. 


3. Nodes of the trial wave function 

The fixed-node approximation is another source of er¬ 
ror within the fixed node DMC framework. During a 
fixed node DMC simulation, the nodes of the trial wave 
function are held fixed to preserve the antisymmetric 
nature of the wave function and avoid collapse to the 
bosonic ground state^^. To accomplish this, in prac¬ 
tice the DMC walkers are prevented from crossing the 
nodes during propagation. If the pseudopotentials are 
sufficiently accurate and other controllable sources of er¬ 
ror such as finite size effects are addressed adequately, 
then the accuracy of the DMC approach will in princi¬ 
ple be limited by the accuracy of the trial wave function 
nodes. If the trial wave function has the exact nodal 
structure, then DMC will project out the exact ground 
state. However, if the nodal structure differs from the 
exact one, then the DMC algorithm will converge to the 
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FIG. 3. DMC energy for ZnSe as a function of the degree of 
exchange a used to generate the trial wave function in DFT. 
For smaller a from 0% < a <~ 50%, the DMC total energies 
exhibit only a small sensitivity. For larger a, the DMC en¬ 
ergy increases by around 0.1 eV/fu at most, indicating that 
the nodal structure is not as good. Total energies are shown 
relative to the minimum calculated energy across all a con¬ 
sidered. 


closest approximation of the ground state subject to the 
constrained nodes. Since the inexact solution has an en¬ 
ergy higher than the exact solution, in principle the nodal 
surface can be optimized by minimizing the total energy 
(although this is challenging in practice). At present, 
a detailed understanding of the magnitude of the fixed 
node error is somewhat lacking, particularly for solids, 
although some recent headway has been made^°. More 
is known about the effects of nodal errors within molec¬ 
ular and/or atomic systems^®’^^ . 

To consider the sensitivity of the final DMC energy on 
the nodal structure. Fig. 3 shows the relative energy of 
ZnSe (eV/fu) obtained from DMC calculations using trial 
wave functions generated from orbitals from hybrid DFT 
calculations, with different degrees of exchange mixing a 
within the PBEIj, framework. The exact nodal structure 
of the trial wave functions is not known, but it is believed 
to span a reasonable range since the orbitals are gener¬ 
ated from theories ranging from DFT-PBE (a = 0%) to 
something similar to Hartree-Fock (a = 100%). For these 
two materials, the main effect of the exchange mixing is 
to modify the relative energies of the p and d orbitals, 
which affects their hybridization and the resultant Kohn- 
Sham orbitals that are used in the Slater determinant. 
Because DMC is variational, the lowest energy obtained 
is the closest upper bound to the exact energy and so we 
minimize with respect to a. This approach of varying 
the exchange weight is often used to obtain better trial 
wave functions for DMC^’^^’^^; some authors vary the N 
parameter in DFT-FC/ to achieve a similar effect^^. 

For the case of ZnSe, from Fig. 3 the sensitivity of 
the total energy to a is very small: for 0% < a < 50% 
the total energies are within error bars and vary by only 
around 0.02 eV/fu. This variation is much smaller than 
the localization errors in the previous section. Only for 
a > 50% does the DMC energy increase, but never by 
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FIG. 4. DMC energy for ZnO as a function of the degree of 
exchange a used to generate the trial wave function in DFT 
using (a) Trail-Needs and (b) Burkatski-Filippi-Dolg pseu¬ 
dopotentials. For both cases there appears to be a minimum 
around a ~ 20% — 30%, indicating a better nodal structure. 
Total energies are shown relative to the minimum calculated 
energy across all a considered. 


more than 0.1 eV/fu for the full range of a considered 
here. Overall, this indicates that for ZnSe the sensitivity 
to the nodal structure of the Slater determinant is small. 
On the other hand for the case of ZnO shown in Fig. 4 
the sensitivity is still small, but for both pseudopoten¬ 
tials there appears to be a minimum somewhere around 
a « 20% — 30%. In some other semiconducting systems 
(MnO^®) we have also found an optimal a « 25% and 
energy variations around 0.1 eV/fu with varying a. This 
is also consistent with other reported findings for transi¬ 
tion metal 3d compound FeO^^. The DMC energies do 
not vary by more than 0.1 eV/fu for 0% < a < 50%, but 
for a = 100% the energy has increased by around 0.3 to 
0.4 eV/fu above the minimum for both pseudopotentials. 
Although the sensitivity of the DMC energy to the nodal 
structure appears to be larger for ZnO (the more ionic 
compound) than ZnSe, it is still small in comparison to 
other effects such as the pseudopotential localization ap¬ 
proximation considered above. 

Regarding nodal errors, it will be interesting in the fu¬ 
ture to establish in detail which solid materials exhibit a 
greater nodal sensitivity and which do not. It will also be 
interesting to assess whether nodal sensitivity is greater 
for excited state wave functions for which the nodal sur¬ 
face may be more complex. Additionally, we emphasize 


that it is still an open question how to most effectively 
further optimize the nodes of the wave function. 


B. Systematically Controllable Approximations 

We now turn to the controllable approximations, for 
which the errors can be made small by converging a simu¬ 
lation parameter. There are several controllable approx¬ 
imations in QMC, which include the finite DMC time 
step, one and many-body finite size effects, the quality 
of the basis set, and the number of configurations. The 
basis set and the number of configurations are relatively 
easy to converge; we consider the more challenging ap¬ 
proximations. We present here results for the DMC time 
step error, the many-body finite size effect for the calcu¬ 
lation of ground state total energies, and the finite size 
effect for the calculation of optical excitation energies. 


1. DMC time step error 


time step (au) 



FIG. 5. Relative DMG energy of ZnSe calculated from an 8 
atom supercell at the F point, as as a function of time step. 
As the time step decreases, the acceptance ratio increases 
and total energy converges towards a fixed value. Stochastic 
uncertainties are smaller than the symbol size. The linear 
extrapolation is shown both for the use of T-moves and the 
localization approximation, demonstrating that the two ap¬ 
proximations give different extrapolated values of the energy. 


In DMC, the time-dependent Schrodinger equation is 
transformed into an imaginary-time diffusion equation 
that includes a source/sink term. A large number of 
“walkers”, which are the possible configurations of all 
N electrons distributed into the 3N dimensional phase 
space (each walker constitutes a 3N dimensional set of 
coordinates) according to the trial wave function. The 
walkers are then propagated according to the dynamics 
of the imaginary time diffusion equation using a Green’s 
function approach. The Green’s function projector is rig¬ 
orously only exact for vanishingly small time step, but 
the propagation of walkers in practice requires a finite 
time step. This finite time step introduces an error in 










7 


the projected energy^^de^ Controlling the time step er¬ 
ror is typically straightforward and can be accomplished 
by performing calculations for a range of time steps and 
extrapolating the result to the limit of zero time step. 
The tradeoff is that smaller time steps require more total 
number of steps for sufficient phase space sampling. If 
the time step is sufficiently small, the results will exhibit 
a linear dependence of the total energy on the time step 
since it the linear term is the first term in the expansion. 

In Fig. 5, the dependence of the total DMC for ZnSe 
on the DMC time step is shown for two cases: when the 
localization approximation is used, and when T-moves 
is used. Both sets of results show the expected trends 
for the dependence of the DMC energy on the timestep; 
we show the linear fit to the calculated results in the 
linear regime. Note that the extrapolated values of the 
DMC energy differ from each other for the two methods. 
In the linear regime the localization approximation ap¬ 
pears to have smaller time step errors than the T-moves 
scheme, although we should note that we did not imple¬ 
ment the recent version of the algorithm^^, which might 
have smaller time step errors. For time steps < 0.01 au, 
the error in the calculated total energy is < 0.1 eV/fu, 
using the algorithms implemented in the QWalk pack¬ 
age. The results for ZnO (not shown) are similar to those 
shown in Fig. 5 for ZnSe. 


(a) 


ZnO, TN pseudopotential 



(b) 


ZnO, BFD pseudopotential 



2. Supercell finite size effects for ground state energies 


FIG. 7. Extrapolation of total energy of ZnO with (a) TN, 
and (b) BFD pseudopotential to thermodynamic limit. Solid 
lines indicate the calculated DMC energies while dashed lines 
include the structure factor correction. 


ZnSe, TN pseudopotential 



FIG. 6. Extrapolation of total energy of ZnSe to thermo¬ 
dynamic limit, using unit cells of different shapes and sizes. 
Solid lines indicate the calculated DMC energies while dashed 
lines include the structure factor correction. 

Simulations of solids often invoke periodic boundary 
conditions applied to a unit cell or a supercell of the 
solid materials. Finite-size effects in QMC can take many 
forms, in particular (i) insufficient twists included for the 
wave function boundary conditions (the analog of insuffi¬ 
cient fc-point sampling in single particle theories) and (ii) 
the many-body finite size effect, which arises from spu¬ 


rious electron-electron image interactions amongst elec¬ 
trons in neighboring cells (this finite size effect is unique 
to many-particle theories and has no analog in single¬ 
particle descriptions). Regarding the former, averaging 
over twists (phases of the wave function) results in a 
faster convergence to the thermodynamic limit as the su¬ 
percell size increases. For all of our supercells we have 
twist-averaged over a 2 x 2 x 2 grid. Since the wave func¬ 
tions for these twists are real, this grid allows for extrap¬ 
olation to larger cells due to the improved computational 
efficiency. 

The many-body finite size effect is of more interest: the 
fictitious electron-electron image interaction introduces 
an artificial correlation that has a stabilizing effect on 
the total energy. This correlation depends on both the 
size and the shape of the supercell and affects the calcu¬ 
lated total energy. Possible correction schemes have been 
discussed in the past^’’“^°, but a systematic investigation 
across a variety of materials carrying out extrapolations 
to large sized supercells is still lacking due to the com¬ 
putational cost of DMC. Here, we illustrate the magni¬ 
tude of the many-body finite size effect, and apply the 
structure-factor S{k) based correction scheme of Chiesa 
et al}^ 

Figs. 6 and 7 show the results for ZnSe and ZnO respec- 











tively. The solid lines are the uncorrected total energies, 
while the dashed lines include the structure factor cor¬ 
rection. In these figures the i-axis shows 1/iV, where N 
is the number of atoms in the supercell, since the lead¬ 
ing order correction for the many-body finite size effect 
scales as 1/iV. In Fig. 6 the blue dots correspond to 
(2 X 2 X I), (2 X 2 X 2), and (2x2x4) supercells of the 
2-atom zincblende fee unit cell. The red dots correspond 
to (2 X 1 X 1), (2 X 2 X I), and (2x2x2) supercells 
built from a 4-atom (wurtzite-like) building block, again 
giving systems of total size 8, 16, and 32 atoms. Finally, 
the green dots correspond to (Ixlxl), (2x1x1), and 
(2x2x1) supercells of the cubic 8 atom conventional 
zincblende unit cell, also giving rise to supercells of 8, 
16, and 32 atoms. In total, for ZnSe we have consid¬ 
ered 9 supercells, 3 each of 8 atoms, 16 atoms, and 32 
atoms. Similarly in Fig. 7 for ZnO the magenta dots 
correspond to (2 x 1 x 1), (2 x 2 x 1), and (2x2x2) 
supercells built from a 4-atom wurtzite building block. 
The blue dots correspond to (1 x 1 x 1), (2 x 1 x 1), and 
(2x2x1) supercells of the 8 atom tetragonal unit cell. 
For both cases, the lines show the best fit linear extrap¬ 
olations of each set of data points to the thermodynamic 
limit TV —>■ oo. The expected 1/TV scaling is observed, 
and all extrapolations yield total ground state energies 
within 0.15 eV/fu of each other. 

According to Figs. 6, 7 the many-body finite size effect 
is large, but the application of the structure factor correc¬ 
tion (dashed lines in Figs. 6, 7) helps substantially. For 
instance for 8 atom supercells of different shapes the total 
energy per fu may be 0.5 -1 eV lower than the extrap¬ 
olated limit for both ZnSe and ZnO. For both however, 
most of the 1 /TV dependence is eliminated and the linear 
fits correspondingly flatten out when the structure factor 
correction is included. In all cases, we find that the use of 
16 or 32 atom supercells, together with S{k) correction, 
are sufficient to resolve energies to within 0.15 eV/fu. 
Further, it is encouraging that all extrapolations, both 
with and without S{k) correction tends towards similar 
values within Ri 0.15 eV/fu of each other. 


3. Supercell finite size effects for optical excitations 

Finally, we consider the calculation of the F-point opti¬ 
cally excited state in ZnSe and ZnO. Since both of these 
semiconductors have a direct gap at F, the calculated 
optical excitation energy will correspond to the optical 
band gap. Accurate calculation of band gaps will be a 
critical piece to establishing the functionality of the QMC 
method, so it is important to assess the scale and mag¬ 
nitude of finite size effects for such systems. To calculate 
the optical excitation energy, we have adopted a proce¬ 
dure that has been used successfully previously^. The 
energy of the first optically excited state is calculated as 

OP = Ev^T-Eg , (3) 



number of atoms 

FIG. 8. Optical excitation energy computed for different su¬ 
percell sizes and shapes for ZnSe. The excitation energy is 
obtained as the difference of the F point ground state and 
first optically excited state energies. 


where Eg is the ground state energy and FTp^r is the 
F-point optically excited state. The energy Er^v is cal¬ 
culated by promoting an electron from the highest occu¬ 
pied Kohn-Sham orbital at F to the lowest unoccupied or¬ 
bital in the construction of the Slater determinant. Here, 
both terms FTp-j-r, Eg are evaluated only at the F point. 
Since the optical excitation energy is an energy differ¬ 
ence between two supercells of the same size and shape, 
the structure factor correction cancels out. Accordingly, 
one might expect that the dominant 1/TV scaling exhib¬ 
ited for the ground state energies in Figs. 6,7 should not 
be present, leaving behind a faster convergence with in¬ 
creasing supercell size. We are unaware of any existing 
detailed analysis of the scaling of finite size effects for 
optical excitation energies. 

To analyze the behavior, in Figs. 8,9 we report the 
computed gaps of ZnSe and ZnO using the same super¬ 
cells that were considered in Figs. 6,7. Since the precise 
scaling of the finite size effect is unknown. Figs. 8,9 show 
the computed gap as a function of supercell size on a lin¬ 
ear/linear scale. The insets show the same data plotted 
on a log/log scale to observe whether the scaling can be 
extracted. For ZnSe (Fig. 8), the TV = 8 atom super¬ 
cells are too small, but it appears that all supercells with 
TV = 16 and TV = 32 are converged within error bars and 
our calculated optical gap of 2.8(2) eV is in agreement 
with experiment. From the inset there is no evidence 
of 1/TV scaling, consistent with the expectation that the 
convergence is faster. 

The situation is more complicated for ZnO (Fig. 9), 
which exhibits the wurtzite structure. It is symmetry- 
broken from cubic zincblende with unequal lattice con¬ 
stants in the c and a, b directions. Correspondingly, there 
are three orbitals at F near the VBM which are slightly 
non-degenerate: one has orbitals aligned along the c axis, 
and two have orbitals aligned within the a, b plane. These 
are plotted in Fig. 10. According to ARPES measure¬ 
ments, the symmetry breaking is small, < 0.1 eV^®. 
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(a) 



(b) 



FIG. 9. Optical excitation energy computed for different su¬ 
percell sizes and shapes for ZnO with both (a) TN and (b) 
BFD pseudopotential. The excitation energy is obtained as 
the difference of the F point ground state and first optically 
excited state energies. The optically excited state trial wave 
functions are constructed by removing an electron both from 
VBM orbitals aligned along the wurtzite a and c axes. 


We address the BFD results for ZnO first, shown in 
Fig. 9b. As for ZnSe the 8 atom supercells are too small, 
but the 16 and 32 atom supercells show more consis¬ 
tent trends, regardless of the supercell used. For the 
32 atom supercells, our calculated excitation energy is 
3.8(2) eV, which slightly overestimates the experimental 
value of 3.4 eV. In addition, the log/log plot in the inset 
also has no evidence of 1/N scaling, suggesting a faster 
convergence. We also tested the sensitivity of our cal¬ 
culated DMC optical excitation energies to the orbital 
from which the electron is removed in the calculation of 
the term Ey^t hr Eq. (3). For the BFD pseudopoten¬ 
tial, we find that (as expected) whether the electron is 
removed from the state with orbitals aligned parallel to 
the c or the a axis, our computed excitation energies are 
not sensitive. Overall, our calculated gaps for ZnO using 
the BFD pseudopotential overestimate the experimental 
value (3.4 eV) by a few tenths of an eV. These results 
are in line with self-consistent GW calculations, which 



FIG. 10. Gharge density of the DFT-PBEO VBM orbitals 
at the F point, using the TN pseudopotential, (a) orbital 
is aligned in the a, b plane, and (b) orbital is parallel to the 
wurtzite c axis. 


find an GW gap of 3.84 eV^®, although we should note 
that non self-consistent GW has a strong sensitivity to 
the calculation starting point for ZnO. 

The situation is different for the TN pseudopotential, 
however, for which our results are shown in Fig. 9a. The 
optical excitation energy has been calculated using the 
same N = 8, 16, 32 supercells that we used for the BFD 
pseudopotential. Compared to the BFD results, the 
trends in the computed TN optical excitations are less 
clear, and it appears that the finite size effects are not 
yet completely converged. Additionally there is an unex¬ 
pected sensitivity of the excitation energy to the orbital 
from which the electron is removed in the calculation of 
Ar-s-r. The reason for the different trends in Fig 9a, b 
is not obvious, but the only difference between the two 
sets results is the pseudopotentials used for the calcula¬ 
tions. We conclude that the Burkatski-Fillipi-Dolg pseu¬ 
dopotential is better able to describe the optically excited 
state wave function. Although there may be other differ¬ 
ences as well, a key difference between the TN and the 
BFD pseudopotentials is the size of the core (Ar core vs. 
Ne core), so it is possible that the ionic cores may be more 
largely perturbed for the excited states. It is interesting 
that throughout this work the differences between the 
two pseudopotentials largely arose in the excited state 
simulations. 


C. Conclusion 

In conclusion, we have carried out a detailed assess¬ 
ment of errors and uncertainties for the simulation of 
ZnSe and ZnO, arising from both controllable and un¬ 
controllable approximations, within the fixed node dif¬ 
fusion Monte Carlo framework. We find that the both 
Trail-Needs and Burkatski-Filippi-Dolg Dirac-Fock pseu¬ 
dopotentials do an excellent job of the calculation of the 
atomization energy. Localization errors introduced by 
non-local pseudopotentials can be very significant, par¬ 
ticularly for pseudopotentials with large cutoff radii. For 
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the ground state the magnitude of the fixed node error, 
as assessed by varying the nodes of the trial wave func¬ 
tion, is comparatively small. Supercell finite size effects 
can have a significant effect on calculated ground state 
energies and extrapolation to the infinite supercell limit 
appears feasible by sampling supercells up to 16 or 32 
atoms in size, particularly with the application of the 
S{k) structure factor correction. When carried out sys¬ 
tematically using the approach outlined here, the result¬ 
ing atomization energies for ZnO and ZnSe are within 0.1 
eV/atom of the experimental value, which is promising 
for the prediction of stability. 

For the calculation of optically excited states, assess¬ 
ment across a variety of supercell sizes and shapes is nec¬ 
essary. We observe that the computed excitations con¬ 
verge quickly with system size, with a scaling that is 
probably faster than linear. We find that the BFD pseu¬ 
dopotential is able to give a good description of the gap of 
ZnO. On the other hand, while the TN pseudopotential 
performed well for the gap of ZnSe, the results are not 
as clear for ZnO. These observations suggest that differ¬ 
ences between pseudopotentials may be more prevalent 
excited state calculations even if not for the ground state, 
and that when transition metals are involved, small core 
pseudopotentials may be important for achieving high 
accuracy. 


The detailed assessment presented here is expected to 
be of use to the QMC modeling community, towards the 
establishment of QMC “best practices” for the simula¬ 
tions of solid materials. We also expect this to provide a 
foundation for further studies that exploit the accuracy 
available in QMC methods. 
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